Progetto

Dataset e Data Cleaning

Il dataset analizzato è “Descrizione tracciato record dati comunali” dell’ISTAT sul numero di decessi in Italia per i mesi da Gennaio ad Aprile e per gli anni dal 2015 al 2020.
La prima modifica effettuata al dataset è stato unpivot dei dati per la classe sesso_anno. Sucessivamente le informazioni riguardo al genere e all’anno di raccolta del dato sono state separate in due differenti colonne, sono state create le colonne DATA, MESE e GIORNO. Infine sono state eliminate le colonne non più necessarie.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)   
library(ggplot2)
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:dplyr':
## 
##     intersect, setdiff, union

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(vctrs)
library(fitdistrplus)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: survival

## Loading required package: npsurv

## Loading required package: lsei

library(actuar)

## 
## Attaching package: 'actuar'

## The following object is masked from 'package:grDevices':
## 
##     cm

df_path <- '/Volumes/HDD_Ale/DSLab/dati-giornalieri-comune/comune_giorno.csv'
data_path_pop_tot <- "/Volumes/HDD_Ale/ds_lab/dataset/tot_pop.csv"
# df_path <- 'C:\\Users\\stefa\\Desktop\\dsLab\\local\\dataset\\comune_giorno.csv'
# data_path_pop_tot <- "C:\\Users\\stefa\\Desktop\\dsLab\\local\\dataset\\tot_pop.csv"
df <-  read.csv(df_path)
df_pop <- read.csv(data_path_pop_tot)
# sistemo il dataset in forma più comoda
df <- df %>%
  gather(key = "SESSO_ANNO", value = "DECESSI", MASCHI_15:FEMMINE_20)
df %>% separate(SESSO_ANNO, c("SESSO", "ANNO"), "_") -> df
df %>% mutate(DATA = as.Date(paste0("0", GE, "20", ANNO), format = "%m%d%Y")) -> df
df$GIORNO <- as.numeric(substr(as.character(df$GE), 2, 3))
df$MESE <- as.numeric(substr(as.character(df$GE), 1, 1))
df$GE <- NULL
df$TOTALE_15 <- NULL
df$TOTALE_16 <- NULL
df$TOTALE_17 <- NULL
df$TOTALE_18 <- NULL
df$TOTALE_19 <- NULL
df$TOTALE_20 <- NULL
summary(df)

Analisi Differenza nelle distribuzioni

L’analisi dei decessi è focalizzata sulla regione Lombardia.
La prima analisi effettuata per dimostrare che c’è stato un incremento significativo dei decessi nel 2020 rispetto agli anni precedenti è un’analisi della varianza tra i gruppi anno 15, 16, 17, 18, 19, 20.
I dati relativi ad Aprile sono stati inoltre eliminati in quanto non totalmente disponibili nell’anno 2020.
La colonna DECESSI riporta diversi valori pari a 9999. Questo dato è un valore dummy per indicare la mancanza del dato in fase di raccolta. Essendo che vogliamo dimostrare la differenza in media delle distribuzioni tra i gruppi consideriamo solamente i comuni che non presentano questa anomalia.
Il dataset così pulito è formato da 622 comuni Lombardi.

df %>% filter(NOME_REGIONE == 'Lombardia',
              DECESSI != 0,
              MESE %in% c(1,2,3))%>% 
  dplyr::select(NOME_COMUNE,ANNO,MESE,DECESSI,DATA) -> df_comuni

#summary(df_comuni)

df_comuni[df_comuni$DECESSI == 9999,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_9999

df_comuni_NO_9999 <- subset(df_comuni, !(df_comuni$NOME_COMUNE %in% comuni_9999$NOME_COMUNE)) 
summary(df_comuni_NO_9999)

##   NOME_COMUNE        ANNO                MESE          DECESSI      
##  Milano : 8522   Length:110978      Min.   :1.000   Min.   : 1.000  
##  Brescia: 2996   Class :character   1st Qu.:1.000   1st Qu.: 1.000  
##  Bergamo: 2196   Mode  :character   Median :2.000   Median : 1.000  
##  Monza  : 1916                      Mean   :2.038   Mean   : 1.188  
##  Como   : 1511                      3rd Qu.:3.000   3rd Qu.: 1.000  
##  Cremona: 1467                      Max.   :3.000   Max.   :17.000  
##  (Other):92370                                                      
##       DATA           
##  Min.   :2015-01-01  
##  1st Qu.:2016-02-29  
##  Median :2018-01-24  
##  Mean   :2017-11-03  
##  3rd Qu.:2019-03-21  
##  Max.   :2020-03-31  
##

#622 comuni senza 9999
length(unique(df_comuni_NO_9999$NOME_COMUNE))

## [1] 622

Aggregando per comune ed anno, sommando il numero di decessi scopriamo la presenza di forti outliers. Questo suggerisce di differenziare l’analisi in categorie differenti, in modo da evitare che le code pesino eccessivamente sulle distribuzioni.

#Tutti mesi (1,2,3)
#61.5 = 30 + (30-9)*1.5
df_aggr <- plyr::ddply(df_comuni_NO_9999, c("NOME_COMUNE","ANNO"), function(x) colSums(x[c("DECESSI")]))
summary(df_aggr)

##                 NOME_COMUNE       ANNO              DECESSI       
##  Acquafredda          :   6   Length:3726        Min.   :   1.00  
##  Acquanegra Cremonese :   6   Class :character   1st Qu.:   9.00  
##  Acquanegra sul Chiese:   6   Mode  :character   Median :  15.00  
##  Adrara San Martino   :   6                      Mean   :  35.39  
##  Agnadello            :   6                      3rd Qu.:  30.00  
##  Agrate Brianza       :   6                      Max.   :4224.00  
##  (Other)              :3690

Per far ciò si considera quindi la classe dei comuni con bassa mortalità, utilizzando come punto di riferimento il massimo $Q_3 + 1.5\cdot(Q_3-Q_1) = 61.5$.

df_aggr[df_aggr$DECESSI > 61,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE),]$NOME_COMUNE -> comuni_bassi
df_aggr_bassi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_bassi)
summary(df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi),])
#188.4 = 98.75 + (98.75-39)*1.5
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi) & df_aggr$DECESSI > 188,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE) & !(df_aggr$NOME_COMUNE %in% comuni_bassi),]$NOME_COMUNE -> comuni_medi
df_aggr_medi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_medi)
df_aggr_alti <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE)

Essendo ancora elevata la presenza di outliers si procede con la classe dei comuni con media mortalità, utilizzando come punto di riferimento il massimo $Q_3 + 1.5(Q_3-Q_1) = 188.4$. I restanti comuni formeranno la classe ad alta mortalità. Infine i comuni di Milano, Bergamo e Brescia verranno considerati a parte in quanto il numero di decessi è molto più elevato della distribuzione dell’ultima classe costruita.

par(mfrow=c(1,1))
boxplot(DECESSI~ANNO,data=df_aggr_bassi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Basso numero di morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")

boxplot(DECESSI~ANNO,data=df_aggr_medi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Medio numero di Morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")

boxplot(DECESSI~ANNO,data=df_aggr_alti[!(df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia")),], outline = TRUE, main="# Decessi nei Comuni Lombardi con Alto numero di Morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")

boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia"),], outline = TRUE, main="# Decessi nel Comune di Milano, Bergamo e Brescia",
        xlab="Anni", ylab="# Decessi", col="lightblue")

Dall’analisi grafica sui boxplot possiamo notare la presenza di un’incremento nell’anno 2020. Per approfondire le cause dell’aumento del numero di decessi si osservano i dati mensili separatamente, applicando la logica delle classi di mortalità.

#Gennaio----------------------------------------------------------------------------------------------------------------------------
df_aggr <- plyr::ddply(df_comuni_NO_9999[df_comuni_NO_9999$MESE == 1,], c("NOME_COMUNE","ANNO"), function(x) colSums(x[c("DECESSI")]))
#summary(df_aggr)
#20.5 = 10 + (10-3)*1.5
df_aggr[df_aggr$DECESSI > 20,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE),]$NOME_COMUNE -> comuni_bassi
df_aggr_bassi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_bassi)
#summary(df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi),])
#70.5 = 39 + (39-18)*1.5
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi) & df_aggr$DECESSI > 70,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE) & !(df_aggr$NOME_COMUNE %in% comuni_bassi),]$NOME_COMUNE -> comuni_medi
df_aggr_medi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_medi)
df_aggr_alti <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE)
par(mfrow=c(1,1))
boxplot(DECESSI~ANNO,data=df_aggr_bassi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Basso numero di morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")

boxplot(DECESSI~ANNO,data=df_aggr_medi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Medio numero di Morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")

boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE != "Milano",], outline = TRUE, main="# Decessi nei Comuni Lombardi con Alto numero di Morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")

boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE == "Milano",], outline = TRUE, main="# Decessi nel Comune di Milano",
        xlab="Anni", ylab="# Decessi", col="lightblue")

Possiamo notare che graficamente non risultano differenze tra le distribuzioni. In particolare risulta un incremento, di cui andrebbe indagata la significatività nel 2017.

#Febbraio----------------------------------------------------------------------------------------------------------------------------
df_aggr <- plyr::ddply(df_comuni_NO_9999[df_comuni_NO_9999$MESE == 2,], c("NOME_COMUNE","ANNO"), function(x) colSums(x[c("DECESSI")]))
#summary(df_aggr)
#18 = 9 + (9-3)*1.5
df_aggr[df_aggr$DECESSI > 18,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE),]$NOME_COMUNE -> comuni_bassi
df_aggr_bassi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_bassi)
#summary(df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi),])
#60 = 33 + (33-15)*1.5
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi) & df_aggr$DECESSI > 60,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE) & !(df_aggr$NOME_COMUNE %in% comuni_bassi),]$NOME_COMUNE -> comuni_medi
df_aggr_medi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_medi)
df_aggr_alti <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE)
par(mfrow=c(2,2))
boxplot(DECESSI~ANNO,data=df_aggr_bassi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Basso numero di morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_medi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Medio numero di Morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE != "Milano",], outline = TRUE, main="# Decessi nei Comuni Lombardi con Alto numero di Morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE == "Milano",], outline = TRUE, main="# Decessi nel Comune di Milano",
        xlab="Anni", ylab="# Decessi", col="lightblue")

Anche per Febbraio non vediamo particolari differenze tra il 2020 e gli altri anni. Ci aspettiamo dunque che tutta la differenza notata nelle distribuzioni totali sia apportata dal mese di Marzo.

#Marzo----------------------------------------------------------------------------------------------------------------------------
df_aggr <- plyr::ddply(df_comuni_NO_9999[df_comuni_NO_9999$MESE == 3,], c("NOME_COMUNE","ANNO"), function(x) colSums(x[c("DECESSI")]))
#summary(df_aggr)
#25,5 = 12 + (12-3)*1.5
df_aggr[df_aggr$DECESSI > 25,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE),]$NOME_COMUNE -> comuni_bassi
df_aggr_bassi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_bassi)
#summary(df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi),])
#69,5 = 32 + (32-7)*1.5
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi) & df_aggr$DECESSI > 69,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE) & !(df_aggr$NOME_COMUNE %in% comuni_bassi),]$NOME_COMUNE -> comuni_medi
df_aggr_medi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_medi)
df_aggr_alti <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE)
par(mfrow=c(2,2))
boxplot(DECESSI~ANNO,data=df_aggr_bassi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Basso numero di morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_medi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Medio numero di Morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[!(df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia")),], outline = TRUE, main="# Decessi nei Comuni Lombardi con Alto numero di Morti",
        xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia"),], outline = TRUE, main="# Decessi nel Comune di Milano, Bergamo e Brescia",
        xlab="Anni", ylab="# Decessi", col="lightblue")

Come intuito precedentemente il mese di Marzo è quello che più influisce sulla differenza tra gli anni nel numero di decessi. Quindi l’analisi si focalizzerà sui dati relativi a Marzo.
Ricordando le ipotesi per il test dell’ANOVA: \[ \begin{itemize} \item[1] I soggetti sono tra loro indipendenti \item[2] Normalità della distribuzione \item[3] Omoschedasticità \end{itemize} \] La 1 è dimostrata per definizione in quanto i soggetti deceduti non verranno ri-analizzati.
L’ipotesi 2 è violata. Già dall’istogramma si nota che le distribuzioni sono non normali.

Per identificare la distribuzione del numero di decessi nei comuni lombardi è stata eseguita una statistica Goodness of Fit con stima della massima verosomiglianza dei parametri. Sono state considerate le distribuzioni teoriche seguenti: \[ \begin{itemize} \item[1] Weibull \item[2] Log Normale \item[3] Gamma \item[1] Log Logistic \item[2] Pareto \end{itemize} \]

#hp:
#   1) soggetti sono tra loro indipendenti. 
#   2) normalità della distribuzione (shapiro.test(x) #pvalue > 0.05)
#   3) omoschedasticità (bartlett.test(dati, gruppi #pvalue>0.5)



#df_aggr_bassi
#   1) Rispettata perchè la persona deceduta non viene più considerata
#   2) Le distribuzioni sono chiaramente NON normali.


gos_stat <- function(x, anno){
  fw <- fitdist(x, "weibull")
  fln <- fitdist(x, "lnorm")
  fg <- fitdist(x, "gamma")
  fll <- fitdist(x, "llogis",start = list(shape = 1, scale = 500))
  fp <- fitdist(x, "pareto",start = list(shape = 1, scale = 500))
  #fb <- fitdist(x, "burr",start = list(shape1 = 0.3, shape2 = 1, rate = 1))
  denscomp(list(fw, fln, fg, fll, fp), main = paste0("Histogram and theoretical densities, year 20",anno, sep=""))
  print(fln)
  return(gofstat(list(fw, fln, fg, fll, fp), fitnames = c("Weibull", "Lognormal", "Gamma", "Loglogistic", "Pareto")))
}

par(mfrow = c(3, 2))
for(i in 15:20) {
  x <- df_aggr_bassi[df_aggr_bassi$ANNO == i,]
  print(paste0("Bassi - Anno 20",i, sep=""))
  print(gos_stat(x$DECESSI,i))
}

par(mfrow = c(3, 2))
for(i in 15:20) {
  x <- df_aggr_medi[df_aggr_medi$ANNO == i,]
  print(paste0("Medi - Anno 20",i, sep=""))
  print(gos_stat(x$DECESSI,i))
}

par(mfrow = c(3, 2))
for(i in 15:20) {
  x <- df_aggr_alti[df_aggr_alti$ANNO == i & !(df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia")),]
  print(paste0("Alti - Anno 20",i, sep=""))
  print(gos_stat(x$DECESSI,i))
}

A titolo esemplificativo prendiamo i comuni con mortalità media e notiamo che la distribuzione Log Normale è quella che mediamente minimizza le statistiche KS, CVM e AD. In particolare, quest’ultima è la più sensibile alle code della distribuzione e la Log Normale risulta quella più stabile anche per questa statistica. Infine notiamo che i coefficienti sitmati sono
$$

par(mfrow = c(3, 2))
for(i in 15:20) {
  x <- df_aggr_medi[df_aggr_medi$ANNO == i,]
  print(paste0("Medi - Anno 20",i, sep=""))
  print(gos_stat(x$DECESSI,i))
}

## [1] "Medi - Anno 2015"

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters:
##          estimate Std. Error
## meanlog 2.2537817 0.06361485
## sdlog   0.7225258 0.04498210
## Goodness-of-fit statistics
##                                Weibull  Lognormal      Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.1065330 0.09847301 0.09556002  0.06917855
## Cramer-von Mises statistic   0.2891741 0.10406118 0.18862769  0.08570355
## Anderson-Darling statistic   1.6267770 0.88629287 1.09581364  0.66380876
##                                 Pareto
## Kolmogorov-Smirnov statistic 0.2319129
## Cramer-von Mises statistic   1.4088299
## Anderson-Darling statistic   7.8402495
## 
## Goodness-of-fit criteria
##                                 Weibull Lognormal    Gamma Loglogistic   Pareto
## Akaike's Information Criterion 867.2280  867.7113 862.7068    865.0353 903.6055
## Bayesian Information Criterion 872.9476  873.4309 868.4264    870.7549 909.3251
## [1] "Medi - Anno 2016"

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters:
##          estimate Std. Error
## meanlog 2.1438531 0.06529743
## sdlog   0.7445053 0.04617188
## Goodness-of-fit statistics
##                                 Weibull  Lognormal      Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.07767132 0.09011304 0.05849273  0.07211654
## Cramer-von Mises statistic   0.10900192 0.14809059 0.06266015  0.09205998
## Anderson-Darling statistic   0.77918414 0.96423074 0.48905673  0.78930128
##                                 Pareto
## Kolmogorov-Smirnov statistic 0.1753318
## Cramer-von Mises statistic   1.1337252
## Anderson-Darling statistic   6.5056750
## 
## Goodness-of-fit criteria
##                                 Weibull Lognormal    Gamma Loglogistic   Pareto
## Akaike's Information Criterion 855.3264  853.6166 850.8753    856.0882 886.5424
## Bayesian Information Criterion 861.0615  859.3517 856.6104    861.8232 892.2775
## [1] "Medi - Anno 2017"

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters:
##          estimate Std. Error
## meanlog 2.1692309 0.06917373
## sdlog   0.7887019 0.04891286
## Goodness-of-fit statistics
##                                Weibull  Lognormal      Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.1059538 0.07482570 0.08889435  0.06783799
## Cramer-von Mises statistic   0.1854201 0.08989504 0.12328535  0.08239378
## Anderson-Darling statistic   1.1421763 0.65550535 0.77660169  0.60159013
##                                 Pareto
## Kolmogorov-Smirnov statistic 0.1822968
## Cramer-von Mises statistic   0.8740724
## Anderson-Darling statistic   5.2279817
## 
## Goodness-of-fit criteria
##                                 Weibull Lognormal    Gamma Loglogistic   Pareto
## Akaike's Information Criterion 876.6234  875.2087 872.6376    876.3367 900.5708
## Bayesian Information Criterion 882.3585  880.9437 878.3727    882.0718 906.3058
## [1] "Medi - Anno 2018"

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters:
##          estimate Std. Error
## meanlog 2.2371787 0.06138730
## sdlog   0.6999229 0.04340698
## Goodness-of-fit statistics
##                                Weibull  Lognormal     Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.1205721 0.07376718 0.1131095  0.07981057
## Cramer-von Mises statistic   0.3592205 0.10822069 0.2780645  0.12439962
## Anderson-Darling statistic   2.0037181 0.63073832 1.4401341  0.73988129
##                                 Pareto
## Kolmogorov-Smirnov statistic 0.2210531
## Cramer-von Mises statistic   1.2986285
## Anderson-Darling statistic   7.5409603
## 
## Goodness-of-fit criteria
##                                 Weibull Lognormal    Gamma Loglogistic   Pareto
## Akaike's Information Criterion 874.1105  861.8263 866.4812    864.5163 906.8896
## Bayesian Information Criterion 879.8455  867.5614 872.2163    870.2514 912.6247
## [1] "Medi - Anno 2019"

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters:
##          estimate Std. Error
## meanlog 2.1756562 0.07311562
## sdlog   0.8304339 0.05170022
## Goodness-of-fit statistics
##                                 Weibull  Lognormal      Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.07941785 0.07714341 0.08283226  0.09222728
## Cramer-von Mises statistic   0.12695871 0.12215267 0.10805279  0.12538874
## Anderson-Darling statistic   0.86175047 0.88274052 0.75445614  0.99702174
##                                 Pareto
## Kolmogorov-Smirnov statistic 0.1602894
## Cramer-von Mises statistic   0.5496101
## Anderson-Darling statistic   3.7537903
## 
## Goodness-of-fit criteria
##                                 Weibull Lognormal    Gamma Loglogistic   Pareto
## Akaike's Information Criterion 881.0664  883.4672 879.3193    889.3455 901.6016
## Bayesian Information Criterion 886.7860  889.1868 885.0389    895.0651 907.3212
## [1] "Medi - Anno 2020"

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters:
##          estimate Std. Error
## meanlog 3.6764143 0.02364755
## sdlog   0.2696236 0.01672031
## Goodness-of-fit statistics
##                                Weibull Lognormal     Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.1418805 0.1113501 0.1260236   0.0934762
## Cramer-von Mises statistic   0.5587452 0.2332009 0.3202473   0.1869218
## Anderson-Darling statistic   3.5173490 1.6167444 2.1178753   1.5650559
##                                  Pareto
## Kolmogorov-Smirnov statistic  0.4699311
## Cramer-von Mises statistic    6.6977298
## Anderson-Darling statistic   32.1513442
## 
## Goodness-of-fit criteria
##                                 Weibull Lognormal    Gamma Loglogistic   Pareto
## Akaike's Information Criterion 1012.928  988.0023 993.1281    996.5303 1229.578
## Bayesian Information Criterion 1018.663  993.7374 998.8632   1002.2653 1235.313

Violata quindi la seconda ipotesi per l’analisi della varianza passiamo a verificare la terza, quella di omeoschedasticità. Notiamo però che il p-value del test di Bartlett è sempre $< 0.05$. Rifiutiamo quindi l’ipotesi nulla di omeoschedasticità.

bartlett.test(df_aggr_bassi$DECESSI, df_aggr_bassi$ANNO)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  df_aggr_bassi$DECESSI and df_aggr_bassi$ANNO
## Bartlett's K-squared = 250.78, df = 5, p-value < 2.2e-16

bartlett.test(df_aggr_medi$DECESSI, df_aggr_medi$ANNO)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  df_aggr_medi$DECESSI and df_aggr_medi$ANNO
## Bartlett's K-squared = 28.332, df = 5, p-value = 3.135e-05

bartlett.test(df_aggr_alti$DECESSI, df_aggr_alti$ANNO)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  df_aggr_alti$DECESSI and df_aggr_alti$ANNO
## Bartlett's K-squared = 19.659, df = 5, p-value = 0.001448

Essendo le due ipotesi per l’ANOVA violate si può procedere con il test non parametrico di Friedman, e con il post-hoc test di Wilcoxon.
Il test ci mostra che c’è almeno una distribuzione differente dalle altre per ogni classe di mortalità, p-value motlo inferiore anche a $0.01$.

## ANOVA per distribuzioni non normali ed eteroschedastiche, con varianze non omogenee,
## Procediamo con il test di Friedman, e con il post-hoc test di Wilcoxon
# https://www.r-statistics.com/2010/02/post-hoc-analysis-for-friedmans-test-r-code/


df_aggr_bassi <- df_aggr_bassi %>% rstatix::convert_as_factor(ANNO)

df_aggr_bassi %>% count(NOME_COMUNE) %>% filter(n < 6) %>% dplyr::select(NOME_COMUNE) -> comuni_del

df_aggr_bassi_Friedman <- df_aggr_bassi[!(df_aggr_bassi$NOME_COMUNE %in% comuni_del$NOME_COMUNE),]

df_aggr_bassi_15 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 15,]$DECESSI
df_aggr_bassi_16 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 16,]$DECESSI
df_aggr_bassi_17 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 17,]$DECESSI
df_aggr_bassi_18 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 18,]$DECESSI
df_aggr_bassi_19 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 19,]$DECESSI
df_aggr_bassi_20 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 20,]$DECESSI


df_aggr_bassi_matrix <- cbind(df_aggr_bassi_15, df_aggr_bassi_16, df_aggr_bassi_17,df_aggr_bassi_18,df_aggr_bassi_19,df_aggr_bassi_20)

#<0.05 differenza significativa tra almeno uno dei gruppi
friedman.test(df_aggr_bassi_matrix)

## 
##  Friedman rank sum test
## 
## data:  df_aggr_bassi_matrix
## Friedman chi-squared = 629.59, df = 5, p-value < 2.2e-16

#Medi
df_aggr_medi <- df_aggr_medi %>% rstatix::convert_as_factor(ANNO)

df_aggr_medi %>% count(NOME_COMUNE) %>% filter(n < 6) %>% dplyr::select(NOME_COMUNE) -> comuni_del

df_aggr_medi_Friedman <- df_aggr_medi[!(df_aggr_medi$NOME_COMUNE %in% comuni_del$NOME_COMUNE),]

df_aggr_medi_15 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 15,]$DECESSI
df_aggr_medi_16 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 16,]$DECESSI
df_aggr_medi_17 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 17,]$DECESSI
df_aggr_medi_18 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 18,]$DECESSI
df_aggr_medi_19 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 19,]$DECESSI
df_aggr_medi_20 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 20,]$DECESSI


df_aggr_medi_matrix <- cbind(df_aggr_medi_15, df_aggr_medi_16, df_aggr_medi_17,df_aggr_medi_18,df_aggr_medi_19,df_aggr_medi_20)

friedman.test(df_aggr_medi_matrix)

## 
##  Friedman rank sum test
## 
## data:  df_aggr_medi_matrix
## Friedman chi-squared = 282.46, df = 5, p-value < 2.2e-16

#Alti
df_aggr_alti <- df_aggr_alti %>% rstatix::convert_as_factor(ANNO)

df_aggr_alti %>% count(NOME_COMUNE) %>% filter(n < 6) %>% dplyr::select(NOME_COMUNE) -> comuni_del

df_aggr_alti_Friedman <- df_aggr_alti[!(df_aggr_alti$NOME_COMUNE %in% comuni_del$NOME_COMUNE),]

df_aggr_alti_15 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 15,]$DECESSI
df_aggr_alti_16 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 16,]$DECESSI
df_aggr_alti_17 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 17,]$DECESSI
df_aggr_alti_18 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 18,]$DECESSI
df_aggr_alti_19 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 19,]$DECESSI
df_aggr_alti_20 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 20,]$DECESSI

block <- 1:length(df_aggr_alti_20)

df_aggr_alti_matrix <- cbind(df_aggr_alti_15, df_aggr_alti_16, df_aggr_alti_17,df_aggr_alti_18,df_aggr_alti_19,df_aggr_alti_20,block)

friedman.test(df_aggr_alti_matrix)

## 
##  Friedman rank sum test
## 
## data:  df_aggr_alti_matrix
## Friedman chi-squared = 106.6, df = 6, p-value < 2.2e-16

Analizzando con il test post-hoc di Wilcox la differenza tra le distribuzioni a due a due notiamo che il $2020$ è statisticamente differente da ogni altro anno.
Abbiamo dei valori significativi anche per altre distribuzioni, come per esempio nei comuni ad alta mortalità il $2017$ è significativamente differente da tutti gli altri anni, come avevamo già notato precedentemente dall’analisi grafica.

#<0.05 differenza significativa

pairwise.wilcox.test(df_aggr_bassi_Friedman$DECESSI,df_aggr_bassi_Friedman$ANNO, paired = TRUE, exact=FALSE, p.adj="none", correct=FALSE)

## 
##  Pairwise comparisons using Wilcoxon signed rank test 
## 
## data:  df_aggr_bassi_Friedman$DECESSI and df_aggr_bassi_Friedman$ANNO 
## 
##    15     16     17     18     19    
## 16 0.020  -      -      -      -     
## 17 0.115  0.798  -      -      -     
## 18 0.692  0.119  0.160  -      -     
## 19 0.782  0.013  0.043  0.685  -     
## 20 <2e-16 <2e-16 <2e-16 <2e-16 <2e-16
## 
## P value adjustment method: none

pairwise.wilcox.test(df_aggr_medi_Friedman$DECESSI,df_aggr_medi_Friedman$ANNO, paired = TRUE, exact=FALSE, p.adj="none", correct=FALSE)

## 
##  Pairwise comparisons using Wilcoxon signed rank test 
## 
## data:  df_aggr_medi_Friedman$DECESSI and df_aggr_medi_Friedman$ANNO 
## 
##    15     16     17     18     19    
## 16 0.048  -      -      -      -     
## 17 0.327  0.263  -      -      -     
## 18 0.636  0.178  0.269  -      -     
## 19 0.434  0.101  0.701  0.743  -     
## 20 <2e-16 <2e-16 <2e-16 <2e-16 <2e-16
## 
## P value adjustment method: none

pairwise.wilcox.test(df_aggr_alti_Friedman$DECESSI,df_aggr_alti_Friedman$ANNO, paired = TRUE, exact=FALSE, p.adj="none", correct=FALSE)

## 
##  Pairwise comparisons using Wilcoxon signed rank test 
## 
## data:  df_aggr_alti_Friedman$DECESSI and df_aggr_alti_Friedman$ANNO 
## 
##    15      16      17      18      19     
## 16 0.6258  -       -       -       -      
## 17 0.0333  0.0357  -       -       -      
## 18 0.4717  0.1580  0.0010  -       -      
## 19 0.1291  0.0591  0.0023  0.7826  -      
## 20 5.2e-08 5.2e-08 5.2e-08 5.3e-08 5.2e-08
## 
## P value adjustment method: none

Possiamo quindi concludere questa parte di analisi confermando che i comuni Lombardi, nel mese di Marzo $2020$, hanno avuto un incremento significativo nel numero dei decessi.

Valutazione dell’incremento

# Prendo la Lombardia
df %>% filter(NOME_REGIONE == 'Lombardia')%>% 
  dplyr::select(NOME_PROVINCIA,NOME_COMUNE,DECESSI,CL_ETA,ANNO,MESE,GIORNO,DATA_INIZIO_DIFF, DATA, SESSO) -> df
#table(df$NOME_COMUNE) # alcuni comuni sono venuti con i nomi sbagliati li sistemo

df$NOME_COMUNE <- as.character(df$NOME_COMUNE)
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Alm\xe8", "Almè")) %>%
  as.data.frame() -> df
# 2
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Barzan\xf2", "Barzanò"))%>%
  as.data.frame() -> df
# 3
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Bascap\xe8", "Bascapè"))%>%
  as.data.frame() -> df
# 7
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Cant\xf9", "Cantù"))%>%
  as.data.frame() -> df
# 11
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Fenegr\xf2", "Fenegrò"))%>%
  as.data.frame() -> df
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Gambol\xf2", "Gambolò"))%>%
  as.data.frame() -> df
#
# 14
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Mont\xf9 Beccaria", "Montù Beccaria"))%>%
  as.data.frame() -> df
# 15
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Muggi\xf2", "Muggiò"))%>%
  as.data.frame() -> df
# 18
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Ro\xe8 Volciano", "Roè Volciano"))%>%
  as.data.frame() -> df
# 19
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Sal\xf2", "Salò"))%>%
  as.data.frame() -> df
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Santa Maria Ho\xe8", "Santa Maria Hoè"))%>%
  as.data.frame() -> df
# 23
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Tem\xf9", "Temù"))%>%
  as.data.frame() -> df
# 25
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Travac\xf2 Siccomario", "Travacò Siccomario"))%>%
  as.data.frame() -> df
# 28
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Vigan\xf2", "Viganò"))%>%
  as.data.frame() -> df
# 29
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Viggi\xf9", "Viggiù"))%>%
  as.data.frame() -> df
# 30
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Villa d'Alm\xe8", "Villa d'Almè"))%>%
  as.data.frame() -> df
# 31
df %>%
  mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Zerbol\xf2", "Zerbolò"))%>%
  as.data.frame() -> df

df_pop <- df_pop[-c(1,3:12,14,15)]
names(df_pop)[2] <- "POPOLAZIONE"
names(df_pop)[1] <- "NOME_COMUNE"

df %>% left_join(df_pop, by='NOME_COMUNE') -> df

## Warning: Column `NOME_COMUNE` joining character vector and factor, coercing into
## character vector

summary(df)
# controllo NA in popolazione
y<-dplyr::filter(df,is.na(POPOLAZIONE))
table(y$NOME_COMUNE)

Decido di sistemare Manualmente per questi comuni la popolazione totale

df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Borgo Mantovano", 5529))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Borgocarbonara", 1931))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Cadrezzate con Osmate", 2659))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Castelgerundo", 1498))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Centro Valle Intelvi", 3227))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Colli Verdi", 1074))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Gornate Olona", 2223))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Piadena Drizzona", 563))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Puegnago del Garda", 3439))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="San Giorgio Bigarello", 13654))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Solbiate con Cagno", 4653))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Torre de' Busi", 2130))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Valvarrone", 551))%>%
  as.data.frame() -> df
df %>%
  mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Vermezzo con Zelo", 5805))%>%
  as.data.frame() -> df

Sono presenti però anche i dati per le provincie(max 3250315 che equivale alla provincia di milano,) decido di sistemare

# Brescia
df %>% filter(POPOLAZIONE!=1265954) ->df
# milano
df %>% filter(POPOLAZIONE!=3250315) -> df
# bergamo
df %>% filter(POPOLAZIONE!=1114590) -> df
# como
df %>% filter(POPOLAZIONE!=599204) -> df
# cremona
df %>% filter(POPOLAZIONE!=358955) -> df
# lecco
df %>% filter(POPOLAZIONE!=337380) -> df
#lodi
df %>% filter(POPOLAZIONE!=230198) -> df
# mantova
df %>% filter(POPOLAZIONE!=412292) -> df
#Pavia
df %>% filter(POPOLAZIONE!=545888) -> df
# sondrio
df %>% filter(POPOLAZIONE!=181095) -> df
# varese
df %>% filter(POPOLAZIONE!=890768) -> df
summary(df)

##                NOME_PROVINCIA   NOME_COMUNE           DECESSI      
##  Milano               :344424   Length:1762884     Min.   :   0.0  
##  Brescia              :240672   Class :character   1st Qu.:   0.0  
##  Bergamo              :238272   Mode  :character   Median :   0.0  
##  Varese               :173100                      Mean   : 712.7  
##  Monza e della Brianza:141540                      3rd Qu.:   0.0  
##  Pavia                :140304                      Max.   :9999.0  
##  (Other)              :484572                                      
##      CL_ETA          ANNO                MESE           GIORNO     
##  Min.   : 0.00   Length:1762884     Min.   :1.000   Min.   : 1.00  
##  1st Qu.:15.00   Class :character   1st Qu.:1.000   1st Qu.: 8.00  
##  Median :17.00   Mode  :character   Median :2.000   Median :16.00  
##  Mean   :16.54                      Mean   :2.415   Mean   :15.55  
##  3rd Qu.:18.00                      3rd Qu.:3.000   3rd Qu.:23.00  
##  Max.   :21.00                      Max.   :4.000   Max.   :31.00  
##                                                                    
##        DATA_INIZIO_DIFF       DATA               SESSO          
##  1 aprile      :906396   Min.   :2015-01-01   Length:1762884    
##  16 aprile     : 99204   1st Qu.:2016-02-28   Class :character  
##  8 aprile      :195672   Median :2017-08-31   Mode  :character  
##  Dati 2020 n.d.:561612   Mean   :2017-08-28                     
##                          3rd Qu.:2019-02-28                     
##                          Max.   :2020-04-30                     
##                          NA's   :3400                           
##   POPOLAZIONE     
##  Min.   :     33  
##  1st Qu.:   4355  
##  Median :   8883  
##  Mean   :  34510  
##  3rd Qu.:  20625  
##  Max.   :1378689  
##

y<-dplyr::filter(df,is.na(DATA))
df %>% anti_join(y, by="DATA") -> df
summary(df)

##                NOME_PROVINCIA   NOME_COMUNE           DECESSI      
##  Milano               :343728   Length:1759484     Min.   :   0.0  
##  Brescia              :240200   Class :character   1st Qu.:   0.0  
##  Bergamo              :237768   Mode  :character   Median :   0.0  
##  Varese               :172812                      Mean   : 714.1  
##  Monza e della Brianza:141260                      3rd Qu.:   0.0  
##  Pavia                :140096                      Max.   :9999.0  
##  (Other)              :483620                                      
##      CL_ETA          ANNO                MESE           GIORNO     
##  Min.   : 0.00   Length:1759484     Min.   :1.000   Min.   : 1.00  
##  1st Qu.:15.00   Class :character   1st Qu.:1.000   1st Qu.: 8.00  
##  Median :17.00   Mode  :character   Median :2.000   Median :16.00  
##  Mean   :16.54                      Mean   :2.416   Mean   :15.52  
##  3rd Qu.:18.00                      3rd Qu.:3.000   3rd Qu.:23.00  
##  Max.   :21.00                      Max.   :4.000   Max.   :31.00  
##                                                                    
##        DATA_INIZIO_DIFF       DATA               SESSO          
##  1 aprile      :904396   Min.   :2015-01-01   Length:1759484    
##  16 aprile     : 99060   1st Qu.:2016-02-28   Class :character  
##  8 aprile      :195192   Median :2017-08-31   Mode  :character  
##  Dati 2020 n.d.:560836   Mean   :2017-08-28                     
##                          3rd Qu.:2019-02-28                     
##                          Max.   :2020-04-30                     
##                                                                 
##   POPOLAZIONE     
##  Min.   :     33  
##  1st Qu.:   4355  
##  Median :   8883  
##  Mean   :  34469  
##  3rd Qu.:  20544  
##  Max.   :1378689  
##

df <- subset(df, !(GIORNO==29 & MESE==2))
df <- subset(df , MESE!=4)

Per comodità eliminiamo tutti

# metto a NA i valori 9999
df[df == 9999] <- NA
df %>% filter(POPOLAZIONE>=60000) -> df_60000 
y<-dplyr::filter(df_60000,is.na(DECESSI))
table(y$NOME_COMUNE)

## 
## Vigevano 
##     1116

# DEVO TOGLIERE VIGEVANO POICHÈ HO 1116 OSSERVAZIONI A NA
df_60000 %>% filter(NOME_COMUNE!='Vigevano') -> df_60000 
table(df_60000$NOME_COMUNE)

## 
##            Bergamo            Brescia      Busto Arsizio  Cinisello Balsamo 
##               9240              10740               7524               7260 
##               Como            Cremona            Legnano             Milano 
##               8268               8064               6348              15552 
##              Monza              Pavia Sesto San Giovanni             Varese 
##               8796               7764               7788               8040

Abbiamno deciso di mostrare le serie storiche di Comuni con più di 60.000 abitanti:

Bergamo
Brescia
1. Busto Arsizio
Cinisello Balsamo
Como
Cremona
Legnano
Milano
Monza
Pavia
Sesto San Giovanni
Varese

Inizialmente prendiamo il dato assoluto dei decessi per cui sostanzialmente andiamo a fare un ranking delle morti con un incremento percentuale maggiore.

listacomuni <- c("Bergamo", "Brescia", "Busto Arsizio", "Cinisello Balsamo", "Como", "Cremona", "Legnano", "Milano", "Monza", "Pavia", "Sesto San Giovanni", "Varese")
datasets <- list()
for(i in listacomuni){
  df_60000 %>% filter(NOME_COMUNE==i) %>%
    mutate(CL_ETA= case_when(CL_ETA<= 6 ~ "GIOVANI",
                             CL_ETA %in% c(7,8,9,10,11,12,13,14)   ~ "MEDI",
                             CL_ETA > 14 ~ "ANZIANI",)) %>% 
    mutate(DECESSI=replace(DECESSI, DECESSI==9999, NA)) -> x
  
  x %>% group_by(DATA) %>%
    summarise(DECESSI = sum(DECESSI))  %>%
    as_tibble() -> x
  assign(paste0('df_comune_',i), x)
  datasets <- append(datasets,paste0('df_comune_',i))
}

# Milano
ts_Bergamo <- zoo(df_comune_Bergamo$DECESSI, df_comune_Bergamo$DATA)
ts_Bergamo <- ts(ts_Bergamo, frequency = 90)
autoplot(ts_Bergamo, main = paste0('Decessi Bergamo 2015-2020'), xlab = 'Anni a 90 giorni', ylab = 'Decessi')

Il codice per il plot si ripete.

ts_total <- ts.union(ts_Bergamo, ts_Brescia, ts_Bs, ts_Cb, ts_Como, ts_Cremona, ts_Legnano, ts_milano, ts_Monza, ts_Pavia, ts_ssg, ts_Varese)
colnames(ts_total)

##  [1] "ts_Bergamo" "ts_Brescia" "ts_Bs"      "ts_Cb"      "ts_Como"   
##  [6] "ts_Cremona" "ts_Legnano" "ts_milano"  "ts_Monza"   "ts_Pavia"  
## [11] "ts_ssg"     "ts_Varese"

colnames(ts_total)[1:12] <- listacomuni
autoplot(ts_total, main = paste0('Decessi Milano'), xlab = 'Marzo', ylab = 'Decessi')

Come possiamo notare Milano in termini assoluti ha avuto l’incremento maggiore per cui nel mese di Marzo ha avuto un rapido incremento fino ad arrivare per vari giorni a 80 decessi al giorno. Da notare anche Bergamo con la media di 30 decessi al giorno nel mese di Marzo.

Decidiamo quindi di creare una variabile che chiameremo “Decessi pro capite” per cui andiamo a stimare le morti giornaliere rispetto alla popolazione (stima del 2019) andando a vedere non più i termini assoluti. Vogliamo vedere se le zone più colpite rimangono le stesse.

listacomuni <- c("Bergamo", "Brescia", "Busto Arsizio", "Cinisello Balsamo", "Como", "Cremona", "Legnano", "Milano", "Monza", "Pavia", "Sesto San Giovanni", "Varese")
datasets <- list()
str(x$DECESSISUPOPOLAZIONE)

## Warning: Unknown or uninitialised column: `DECESSISUPOPOLAZIONE`.

##  NULL

for(i in listacomuni){
  df_60000 %>% filter(NOME_COMUNE==i) %>%
          mutate(CL_ETA= case_when(CL_ETA<= 6 ~ "GIOVANI",
                                   CL_ETA %in% c(7,8,9,10,11,12,13,14)   ~ "MEDI",
                                  CL_ETA > 14 ~ "ANZIANI",)) %>% 
          mutate(DECESSISUPOPOLAZIONE = (DECESSI/POPOLAZIONE)) %>% 
          mutate(DECESSI=replace(DECESSI, DECESSI==9999, NA)) -> x
  
  x %>% group_by(DATA) %>%
    summarise(DECESSISUPOPOLAZIONE = sum(DECESSISUPOPOLAZIONE))  %>%
    as_tibble() -> x
  assign(paste0('df_comune_',i), x)
  datasets <- append(datasets,paste0('df_comune_',i))
}

# Milano
ts_Bergamo <- zoo(df_comune_Bergamo$DECESSISUPOPOLAZIONE, df_comune_Bergamo$DATA)
ts_Bergamo <- ts(ts_Bergamo, frequency = 90)
autoplot(ts_Bergamo, main = paste0('Decessi Bergamo 2015-2020'), xlab = 'Anni a 90 giorni', ylab = 'Decessi')

Il codice per il plot si ripete.

ts_total <- ts.union(ts_Bergamo, ts_Brescia, ts_Bs, ts_Cb, ts_Como, ts_Cremona, ts_Legnano, ts_milano, ts_Monza, ts_Pavia, ts_ssg, ts_Varese)
# colnames(ts_total)
colnames(ts_total)[1:12] <- listacomuni
autoplot(ts_total, main = paste0('Decessi Milano'), xlab = 'Marzo', ylab = 'Decessi')

Al contrario del grafico precedente se consideriamo i morti rispetto alla popolazione notiamo che non è Milano la città con il record più negativo, ma diventa Bergamo che è seguito dal Cremona e Brescia.

Decidiamo di analizzare i NULL, per cui secondo le specifiche del dataset questi comuni non hanno avuto un incremento del 20% della mortalità dell’anno 2020. Perciò decidiamo arbitrariamente di sostituirli con i valori delle medie dei 5 anni precedenti. Così facendo notiamo fin da subito un incremento molto evidente della mortalità per il mese di Marzo anno 2020. Supponiamo perciò si tratti non solo unicamente del Covid poichè è per noio impossibile da dimostrare, ma sicuramente c’è stato un forte impatto.

df %>% filter(ANNO %in% c(15,16,17,18,19),
              MESE %in% c(1,2,3),
              DECESSI != 9999) %>% 
        group_by(NOME_COMUNE,ANNO,MESE) %>% 
        summarise(DECESSI = sum(DECESSI)) %>% group_by(NOME_COMUNE,MESE) %>% 
                  summarise(DECESSI=mean(DECESSI)) -> mean_dataset
df %>% filter(DECESSI == 9999,
              MESE %in% c(1,2,3)) %>% 
        mutate(DECESSI = NA) %>% 
        group_by(NOME_COMUNE,ANNO,MESE) %>% 
        summarise(DECESSI = sum()) %>% 
        left_join(mean_dataset, by = c("NOME_COMUNE" = "NOME_COMUNE", "MESE" = "MESE")) %>% mutate(DECESSI.x = as.integer(DECESSI.y)) -> finale
df  %>%filter(MESE %in% c(1,2,3)) %>% mutate(DECESSI = ifelse(DECESSI==9999, NA, DECESSI)) %>% group_by(NOME_COMUNE,ANNO,MESE) %>% 
        summarise(DECESSI = sum(DECESSI)) -> oprpr
finale_sistemati_null <- oprpr %>% left_join(finale, by = c("NOME_COMUNE" = "NOME_COMUNE", "MESE" = "MESE")) %>% mutate(DECESSI = ifelse(is.na(DECESSI), as.integer(DECESSI.y), DECESSI))
colnames(finale_sistemati_null)[2] <- "ANNO" 
finale_sistemati_null <- finale_sistemati_null[-c(5:7)]
finale_sistemati_null %>%
  group_by(ANNO, MESE) %>% 
  summarise(DECESSI = sum(DECESSI)) -> finale_sistemati_null

Progetto

Alessandro Pontini, Stefano Biondi

5/5/2020

Dataset e Data Cleaning

Analisi Differenza nelle distribuzioni

Valutazione dell’incremento